
###############################################################################
# Senators at Home 
# By Jaclyn Kaslovsky  
###############################################################################

# Load the packages
require(ggplot2)
require(interflex)
require(readstata13)
require(lfe)
require(plyr)
require(dplyr)
require(sf)
require(usmap)
require(ggridges)
require(lemon)
require(stringr)
require(ggpubr)
library(tigris)

###############################
# Figures 1 and 3: Large Maps
###############################

#Load in the data
datamap <- read.dta13("resources_tocounties_final.dta")

#Load in census shapefile
counties_2010 <- counties(cb = TRUE, year=2010) %>%
  shift_geometry()

states_2010 <- states(cb = TRUE, year=2010) %>%
  shift_geometry()


#Get fips codes and remove puerto rico
counties_2010$GEOID <- str_sub(counties_2010$GEO_ID, start= -5)
counties_2010 <- counties_2010[counties_2010$STATE!=72,]
states_2010 <- states_2010[states_2010$STATE!=72,]

#A function to remove map axes
ditch_the_axes <- theme(
  axis.text = element_blank(),
  axis.line = element_blank(),
  axis.ticks = element_blank(),
  panel.border = element_blank(),
  panel.grid = element_blank(),
  axis.title = element_blank()
)


#Load the data to make a map of trips home and collapse it to the county level
tomap <- ddply(datamap, .(fips, stateabrev), summarize, Sum=sum(traveled))
tomap$logtraveled <- log(tomap$Sum +1)
colnames(tomap)[1] <- "GEOID"
tomap_final <- left_join(counties_2010,tomap, by="GEOID")

#Map it 
map_traveled <- tomap_final %>% 
  ggplot() +
  geom_sf(aes(fill=logtraveled), lwd = 0) +
  theme_bw() + ditch_the_axes  + 
  labs(fill = "ln(Number of Trips to County + 1)") +
  scale_fill_gradient(low="grey95", high="midnightblue") +
  theme(legend.position="bottom") + geom_sf(data=states_2010, fill=NA, color="grey60") +
  coord_sf(crs = st_crs(2163), xlim = c(-2500000, 2500000)) 
ggsave("Figure1.jpeg",map_traveled,width = 12, height = 8, dpi=600) 


#Do the same for staff
tomap <- ddply(datamap, .(fips, stateabrev), summarize, Max=max(number_offices, na.rm=TRUE))
tomap$Max[tomap$Max>=1] <- 1
colnames(tomap)[1] <- "GEOID"
tomap_final <- left_join(counties_2010,tomap, by="GEOID")
tomap_final$Max[tomap_final$GEOID=="11001"] <- 0
tomap_final$Max <- as.factor(tomap_final$Max)

map_staffed <- tomap_final %>% 
  ggplot() +
  geom_sf(aes(fill=Max), lwd = 0.5, color="grey80") +
  theme_bw() + ditch_the_axes  + 
  labs(fill = "Counties with At Least One District Office") +   
  scale_fill_manual(values=c("white","midnightblue")) +
  theme(legend.position="bottom") + geom_sf(data=states_2010, fill=NA, color="grey60") +
  coord_sf(crs = st_crs(2163), xlim = c(-2500000, 2500000)) 
ggsave("Figure3.jpeg",map_staffed,width = 12, height = 8, dpi=600) 


########################
# Two Senator Examples #
#########################

###### Visit graphs

#Limit the data to Ron Johnson by congress
datamap_ron112 <- datamap[datamap$SenatorName=="JOHNSON, Ron" & datamap$cong==112,]

tomap <- ddply(datamap_ron112, .(fips), summarize, Sum=sum(traveled))
tomap$logtraveled <- log(tomap$Sum +1)
df <- us_map(regions = "counties")
tomap_final <- left_join(df, tomap, by = "fips") 
all <- tomap_final[,c("fips","logtraveled")]

map112 <- plot_usmap(regions="counties", include = c("WI"), data=all, values="logtraveled", color="grey80")  +
  theme_bw() + ditch_the_axes  + 
  labs(fill = "ln(Number of Trips to County + 1)") +
  scale_fill_gradient(low="grey95", high="red4") +
  theme(legend.position="bottom") + ggtitle("112th Congress")

datamap_ron113 <- datamap[datamap$SenatorName=="JOHNSON, Ron" & datamap$cong==113,]

tomap <- ddply(datamap_ron113, .(fips), summarize, Sum=sum(traveled))
tomap$logtraveled <- log(tomap$Sum +1)
df <- us_map(regions = "counties")
tomap_final <- left_join(df, tomap, by = "fips") 
all <- tomap_final[,c("fips","logtraveled")]

map113 <- plot_usmap(regions="counties", include = c("WI"), data=all, values="logtraveled", color="grey80")  +
  theme_bw() + ditch_the_axes  + 
  labs(fill = "ln(Number of Trips to County + 1)") +
  scale_fill_gradient(low="grey95", high="red4") +
  theme(legend.position="bottom") + ggtitle("113th Congress")

datamap_ron114 <- datamap[datamap$SenatorName=="JOHNSON, Ron" & datamap$cong==114,]

tomap <- ddply(datamap_ron114, .(fips), summarize, Sum=sum(traveled))
tomap$logtraveled <- log(tomap$Sum +1)
df <- us_map(regions = "counties")
tomap_final <- left_join(df, tomap, by = "fips") 
all <- tomap_final[,c("fips","logtraveled")]

map114 <-  plot_usmap(regions="counties", include = c("WI"), data=all, values="logtraveled", color="grey80")  +
  theme_bw() + ditch_the_axes  + 
  labs(fill = "ln(Number of Trips to County + 1)") +
  scale_fill_gradient(low="grey95", high="red4") +
  theme(legend.position="bottom") + ggtitle("114th Congress (election year)")

datamap_ron115 <- datamap[datamap$SenatorName=="JOHNSON, Ron" & datamap$cong==115,]

tomap <- ddply(datamap_ron115, .(fips), summarize, Sum=sum(traveled))
tomap$logtraveled <- log(tomap$Sum +1)
df <- us_map(regions = "counties")
tomap_final <- left_join(df, tomap, by = "fips") 
all <- tomap_final[,c("fips","logtraveled")]

map115 <- plot_usmap(regions="counties", include = c("WI"), data=all, values="logtraveled", color="grey80")  +
  theme_bw() + ditch_the_axes  + 
  labs(fill = "ln(Number of Trips to County + 1)") +
  scale_fill_gradient(low="grey95", high="red4") +
  theme(legend.position="bottom") + ggtitle("115th Congress")

#Combine the maps and save
ron <- ggarrange(map112, map113, map114, map115, ncol=4,  common.legend = TRUE, legend="bottom")
ggsave("Figure2b.pdf",ron,width = 12, height = 6, dpi=600) 


#Now do the same thing for Tammy Baldwin 
datamap_tammy113 <- datamap[datamap$SenatorName=="BALDWIN, Tammy" & datamap$cong==113,]

tomap <- ddply(datamap_tammy113, .(fips), summarize, Sum=sum(traveled))
tomap$logtraveled <- log(tomap$Sum +1)
df <- us_map(regions = "counties")
tomap_final <- left_join(df, tomap, by = "fips") 
all <- tomap_final[,c("fips","logtraveled")]

map113 <- plot_usmap(regions="counties", include = c("WI"), data=all, values="logtraveled", color="grey80")  +
  theme_bw() + ditch_the_axes  + 
  labs(fill = "ln(Number of Trips to County + 1)") +
  scale_fill_gradient(low="grey95", high="dodgerblue4") +
  theme(legend.position="bottom") + ggtitle("113th Congress")  + theme(plot.title = element_text(size=16)) + theme(legend.title=element_text(size=14), legend.text=element_text(size=10))

datamap_tammy114 <- datamap[datamap$SenatorName=="BALDWIN, Tammy" & datamap$cong==114,]

tomap <- ddply(datamap_tammy114, .(fips), summarize, Sum=sum(traveled))
tomap$logtraveled <- log(tomap$Sum +1)
df <- us_map(regions = "counties")
tomap_final <- left_join(df, tomap, by = "fips") 
all <- tomap_final[,c("fips","logtraveled")]

map114 <-  plot_usmap(regions="counties", include = c("WI"), data=all, values="logtraveled", color="grey80")  +
  theme_bw() + ditch_the_axes  + 
  labs(fill = "ln(Number of Trips to County + 1)") +
  scale_fill_gradient(low="grey95", high="dodgerblue4") +
  theme(legend.position="bottom") + ggtitle("114th Congress")  + theme(plot.title = element_text(size=16)) + theme(legend.title=element_text(size=14), legend.text=element_text(size=10))

datamap_tammy115 <- datamap[datamap$SenatorName=="BALDWIN, Tammy" & datamap$cong==115,]

tomap <- ddply(datamap_tammy115, .(fips), summarize, Sum=sum(traveled))
tomap$logtraveled <- log(tomap$Sum +1)
df <- us_map(regions = "counties")
tomap_final <- left_join(df, tomap, by = "fips") 
all <- tomap_final[,c("fips","logtraveled")]

map115 <- plot_usmap(regions="counties", include = c("WI"), data=all, values="logtraveled", color="grey80")  +
  theme_bw() + ditch_the_axes  + 
  labs(fill = "ln(Number of Trips to County + 1)") +
  scale_fill_gradient(low="grey95", high="dodgerblue4") +
  theme(legend.position="bottom") + ggtitle("115th Congress (election year)") + theme(plot.title = element_text(size=16)) + theme(legend.title=element_text(size=14), legend.text=element_text(size=10))

#Combine the graphs and save
tammy <- ggarrange(map113, map114, map115, ncol=3,  common.legend = TRUE, legend="bottom")
ggsave("Figure2a.pdf", tammy,width = 12, height = 6, dpi=600) 

###### Staffing Graphs 

#Ron Johnson 
staff_ron <- datamap[datamap$SenatorName=="JOHNSON, Ron",]

tomap <- ddply(staff_ron, .(fips), summarize, meanstaff=mean(staff_per))
df <- us_map(regions = "counties")
tomap_final <- left_join(df, tomap, by = "fips") 
all <- tomap_final[,c("fips","meanstaff")]
all$max <- all$meanstaff > 0

ronstaffmax <-  plot_usmap(regions="counties", include = c("WI"), data=all, values="max", color="grey80")  +
  theme_bw() + ditch_the_axes  + 
  labs(fill = "") +   
  theme(legend.position="bottom") + ggtitle("Senator Ron Johnson (R-WI)") +
  theme(plot.title = element_text(hjust = 0.5))  + scale_fill_manual(labels = c("No Office", "Office"),values=c("grey95","red4"))

#Tammy Baldwin
staff_tammy <- datamap[datamap$SenatorName=="BALDWIN, Tammy",]
tomap <- ddply(staff_tammy, .(fips), summarize, meanstaff=mean(staff_per))
df <- us_map(regions = "counties")
tomap_final <- left_join(df, tomap, by = "fips") 
all <- tomap_final[,c("fips","meanstaff")]
all$max <- all$meanstaff > 0

tammystaffmax <-  plot_usmap(regions="counties", include = c("WI"), data=all, values="max", color="grey80")  +
  theme_bw() + ditch_the_axes  + 
  labs(fill = "") +   
  theme(legend.position="bottom") + ggtitle("Senator Tammy Baldwin (D-WI)") +
  theme(plot.title = element_text(hjust = 0.5)) + scale_fill_manual(labels = c("No Office", "Office"),values=c("grey95","dodgerblue4"))

#Combine the graphs and save
figure4 <- ggarrange(tammystaffmax, ronstaffmax, ncol=2,nrow=1)
ggsave("Figure4.pdf",figure4,width = 8, height = 8, dpi=600) 


####################################
# Policy Disagreement over Time   #
####################################

#Load the data and create a ridge plot
data <- read.dta13("complete_CCES_final.dta",convert.factors=FALSE)
data$race <- as.factor(data$race)
data <- data[!is.na(data$approval_binary),]
data <- data[!is.na(data$votematchmean2),]
data <- data[!is.na(data$copartisan),]

ridge <- ggplot(data, aes(x = votematchmean2, y = as.factor(year), weight = weight_cumulative)) + 
  geom_density_ridges() + theme_minimal()+ xlab("Policy Disagreement") +
  ylab("Year") 

ggsave("FigureD1.pdf",ridge,width = 8, height = 8, dpi=600) 

#######################################
# Interflex Graphs: Approval: Visits #
######################################

#Create the binary news interest variable
data$news_watch[data$newsint==1 | data$newsint==2] <-1
data$news_watch[data$newsint==3 | data$newsint==4 | data$newsint==7] <- 0
data_msa <- data[!is.na(data$traveled_msa),]

#Subset the data into news watchers and non-news watchers
news <- data[!is.na(data$news_watch),]
news_msa <- data_msa[!is.na(data_msa$news_watch),]

news_vote <- news[!is.na(news$voted_for_inc),]
news_msa_vote <- news_msa[!is.na(news_msa$voted_for_inc),]

#County
visits_county_news <- interflex(estimator="binning",full.moderate = TRUE, Y = "approval_binary", D = "log_traveled", X = "votematchmean2", 
                                Z = c("copartisan", "electionyear", "seniority", "nokken_poole_dim1", "chair", "power", "race", "female", "age"),
                                data = news[news$news_watch==1,],  main = "", weights="weight_cumulative", ylim=c(-0.03, 0.02),
                                FE = c("member_county", "party_year"),vcov.type="cluster", treat.type="continuous", cl="icpsr", na.rm=TRUE, bin.labs=FALSE, Xlabel="Policy Disagreement", 
                                Dlabel="ln(Visits + 1)", Ylabel="Approval", theme.bw=TRUE, Xdistr="density")


visits_county_nonews <- interflex(estimator="binning",full.moderate = TRUE, Y = "approval_binary", D = "log_traveled", X = "votematchmean2", 
                                  Z = c("copartisan", "electionyear", "seniority", "nokken_poole_dim1", "chair", "power", "race", "female", "age"),
                                  data = news[news$news_watch==0,],  main = "", weights="weight_cumulative", ylim=c(-0.03, 0.02),
                                  FE = c("member_county", "party_year"),vcov.type="cluster", treat.type="continuous", cl="icpsr", na.rm=TRUE, bin.labs=FALSE, Xlabel="Policy Disagreement", 
                                  Dlabel="ln(Visits + 1)", Ylabel="Approval", theme.bw=TRUE, Xdistr="density")
#MSA
visits_msa_news <- interflex(estimator="binning",full.moderate = TRUE, Y = "approval_binary", D = "log_traveled_msa", X = "votematchmean2", 
                             Z = c("copartisan", "electionyear", "seniority", "nokken_poole_dim1", "chair", "power", "race", "female", "age"),
                             data = news_msa[news_msa$news_watch==1,],  main = "", weights="weight_cumulative", ylim=c(-0.03, 0.02),
                             FE = c("member_msa", "party_year"), vcov.type="cluster", treat.type="continuous", cl="icpsr", na.rm=TRUE, bin.labs=FALSE, Xlabel="Policy Disagreement", 
                             Dlabel="ln(Visits + 1)", Ylabel="Approval", theme.bw=TRUE, Xdistr="density")


visits_msa_nonews <-interflex(estimator="binning",full.moderate = TRUE, Y = "approval_binary", D = "log_traveled_msa", X = "votematchmean2", 
                              Z = c("copartisan", "electionyear", "seniority", "nokken_poole_dim1", "chair", "power", "race", "female", "age"),
                              data = news_msa[news_msa$news_watch==0,],  main = "", weights="weight_cumulative", ylim=c(-0.03, 0.02),
                              FE = c("member_msa", "party_year"),vcov.type="cluster", treat.type="continuous", cl="icpsr", na.rm=TRUE, bin.labs=FALSE, Xlabel="Policy Disagreement", 
                              Dlabel="ln(Visits + 1)", Ylabel="Approval", theme.bw=TRUE, Xdistr="density")


#Pull out the necessary data from the interflex models to later graph
county_news <- data.frame(visits_county_news$est.bin$`log_traveled=0 (50%)`[,1:5])
county_news$type <- "County"


county_news_linear <- data.frame(visits_county_news$est.lin$`log_traveled=0 (50%)`[,1:5])
county_news_linear$type <- "County"

county_nonews <- data.frame(visits_county_nonews$est.bin$`log_traveled=0 (50%)`[,1:5])
county_nonews$type <- "County"

county_nonews_linear <- data.frame(visits_county_nonews$est.lin$`log_traveled=0 (50%)`[,1:5])
county_nonews_linear$type <- "County"


msa_news <- data.frame(visits_msa_news$est.bin$`log_traveled_msa=0.693 (50%)`[,1:5])
msa_news$type <- "Metropolitan Statistical Area"
msa_news_linear <- data.frame(visits_msa_news$est.lin$`log_traveled_msa=0.693 (50%)`[,1:5])
msa_news_linear$type <- "Metropolitan Statistical Area"

msa_nonews <- data.frame(visits_msa_nonews$est.bin$`log_traveled_msa=0.693 (50%)`[,1:5])
msa_nonews$type <- "Metropolitan Statistical Area"
msa_nonews_linear <- data.frame(visits_msa_nonews$est.lin$`log_traveled_msa=0.693 (50%)`[,1:5])
msa_nonews_linear$type <- "Metropolitan Statistical Area"

#Combine the dataframes
total <- bind_rows(county_news, msa_news)
total_nonews <- bind_rows(county_nonews, msa_nonews)

total_linear <- bind_rows(county_news_linear, msa_news_linear)
total_linear_nonews <- bind_rows(county_nonews_linear, msa_nonews_linear)


#Graph it
visits_approval_news <- ggplot(total_linear, aes(X)) +
  geom_line(aes(y = ME, color=factor(type)), size=0.8,alpha = .12, position=position_dodge(width=.1), linetype=1) +
  geom_ribbon(aes(ymin = lower.CI.95.., ymax =upper.CI.95..,fill=type), alpha = .12, position=position_dodge(width=.1)) +
  geom_errorbar(data=total,size=1, aes(x=x0,ymin=CI.lower,ymax=CI.upper,color=type),
                width=0,position=position_dodge(width=.1))+
  geom_point(data=total,aes(x=x0,y=coef,fill=type,color=type), size=2, position=position_dodge(width=.1))  +
  geom_hline(yintercept=0,linetype=2, size=.8) +
  scale_color_manual(values=c("plum4", "black"), name="") +
  scale_fill_manual(values=c("plum4", "black"), name="") + ylim(-.03,.02) +
  ylab("Marginal Effect of ln(Visits + 1) on Approval")+ xlab("Policy Disagreement") + theme_minimal() +
  theme(legend.position="bottom",
        plot.title = element_text(hjust = 0.5))  + ggtitle("Regular News Followers") 


visits_approval_nonews <- ggplot(total_linear_nonews, aes(X)) +
  geom_line(aes(y = ME, color=factor(type)), size=0.8,alpha = .12, position=position_dodge(width=.1), linetype=1) +
  geom_ribbon(aes(ymin = lower.CI.95.., ymax =upper.CI.95..,fill=type), alpha = .12, position=position_dodge(width=.1)) +
  geom_errorbar(data=total_nonews,size=1, aes(x=x0,ymin=CI.lower,ymax=CI.upper,color=type),
                width=0,position=position_dodge(width=.1))+
  geom_point(data=total_nonews,aes(x=x0,y=coef,fill=type,color=type), size=2, position=position_dodge(width=.1))  +
  geom_hline(yintercept=0,linetype=2, size=.8) +
  scale_color_manual(values=c("plum4", "black"), name="") +
  scale_fill_manual(values=c("plum4", "black"), name="") +ylim(-.03,.02) +
  ylab("Marginal Effect of ln(Visits + 1) on Approval")+ xlab("Policy Disagreement") + theme_minimal() +
  theme(legend.position="none")+
  theme(axis.text.y = element_blank(),
        axis.title.y= element_blank(),
        plot.title = element_text(hjust = 0.5)) + ggtitle("Non-Regular News Followers") 

visits_approval <- grid_arrange_shared_legend(visits_approval_news, visits_approval_nonews,  ncol=2)
ggsave("Figure5.pdf", plot=visits_approval,width = 12, height = 5, dpi=600)

####################################
#Interflex Graphs: Approval: Staff # 
####################################

data2 <- data[!is.na(data$staff_per),]
data2_msa <- data2[!is.na(data2$staff_per_msa),]
news <- data2[!is.na(data2$news_watch),]
news_msa <- data2_msa[!is.na(data2_msa$news_watch),]

#County
staff_county_news <- interflex(estimator="binning",full.moderate = TRUE, Y = "approval_binary", D = "log_staffed", X = "votematchmean2", 
                                Z = c("copartisan", "electionyear", "seniority", "nokken_poole_dim1", "chair", "power", "race", "female", "age"),
                                data = news[news$news_watch==1,],  main = "", weights="weight_cumulative",
                                FE = c("member_county", "party_year"),vcov.type="cluster", treat.type="continuous", cl="icpsr", na.rm=TRUE, bin.labs=FALSE, Xlabel="Policy Disagreement", 
                                Dlabel="ln(staff + 1)", Ylabel="Approval", theme.bw=TRUE, Xdistr="density")


staff_county_nonews <- interflex(estimator="binning",full.moderate = TRUE, Y = "approval_binary", D = "log_staffed", X = "votematchmean2", 
                                  Z = c("copartisan", "electionyear", "seniority", "nokken_poole_dim1", "chair", "power", "race", "female", "age"),
                                  data = news[news$news_watch==0,],  main = "", weights="weight_cumulative",
                                  FE = c("member_county", "party_year"),vcov.type="cluster", treat.type="continuous", cl="icpsr", na.rm=TRUE, bin.labs=FALSE, Xlabel="Policy Disagreement", 
                                  Dlabel="ln(staff + 1)", Ylabel="Approval", theme.bw=TRUE, Xdistr="density")
#MSA
staff_msa_news <- interflex(estimator="binning",full.moderate = TRUE, Y = "approval_binary", D = "log_staffed_msa", X = "votematchmean2", 
                             Z = c("copartisan", "electionyear", "seniority", "nokken_poole_dim1", "chair", "power", "race", "female", "age"),
                             data = news_msa[news_msa$news_watch==1,],  main = "", weights="weight_cumulative", 
                             FE = c("member_msa", "party_year"), vcov.type="cluster", treat.type="continuous", cl="icpsr", na.rm=TRUE, bin.labs=FALSE, Xlabel="Policy Disagreement", 
                             Dlabel="ln(staff + 1)", Ylabel="Approval", theme.bw=TRUE, Xdistr="density")


staff_msa_nonews <-interflex(estimator="binning",full.moderate = TRUE, Y = "approval_binary", D = "log_staffed_msa", X = "votematchmean2", 
                              Z = c("copartisan", "electionyear", "seniority", "nokken_poole_dim1", "chair", "power", "race", "female", "age"),
                              data = news_msa[news_msa$news_watch==0,],  main = "", weights="weight_cumulative", 
                              FE = c("member_msa", "party_year"),vcov.type="cluster", treat.type="continuous", cl="icpsr", na.rm=TRUE, bin.labs=FALSE, Xlabel="Policy Disagreement", 
                              Dlabel="ln(staff + 1)", Ylabel="Approval", theme.bw=TRUE, Xdistr="density")


#Pull out the necessary data from the interflex models to later graph
county_news <- data.frame(staff_county_news$est.bin$`log_staffed=0 (50%)`[,1:5])
county_news$type <- "County"


county_news_linear <- data.frame(staff_county_news$est.lin$`log_staffed=0 (50%)`[,1:5])
county_news_linear$type <- "County"

county_nonews <- data.frame(staff_county_nonews$est.bin$`log_staffed=0 (50%)`[,1:5])
county_nonews$type <- "County"

county_nonews_linear <- data.frame(staff_county_nonews$est.lin$`log_staffed=0 (50%)`[,1:5])
county_nonews_linear$type <- "County"


msa_news <- data.frame(staff_msa_news$est.bin$`log_staffed_msa=2.285 (50%)`[,1:5])
msa_news$type <- "Metropolitan Statistical Area"
msa_news_linear <- data.frame(staff_msa_news$est.lin$`log_staffed_msa=2.285 (50%)`[,1:5])
msa_news_linear$type <- "Metropolitan Statistical Area"

msa_nonews <- data.frame(staff_msa_nonews$est.bin$`log_staffed_msa=2.246 (50%)`[,1:5])
msa_nonews$type <- "Metropolitan Statistical Area"
msa_nonews_linear <- data.frame(staff_msa_nonews$est.lin$`log_staffed_msa=2.246 (50%)`[,1:5])
msa_nonews_linear$type <- "Metropolitan Statistical Area"

#Combine the dataframes
total <- bind_rows(county_news, msa_news)
total_nonews <- bind_rows(county_nonews, msa_nonews)

total_linear <- bind_rows(county_news_linear, msa_news_linear)
total_linear_nonews <- bind_rows(county_nonews_linear, msa_nonews_linear)


#Graph it
staff_approval_news <- ggplot(total_linear, aes(X)) +
  geom_line(aes(y = ME, color=factor(type)), size=0.8,alpha = .12, position=position_dodge(width=.1), linetype=1) +
  geom_ribbon(aes(ymin = lower.CI.95.., ymax =upper.CI.95..,fill=type), alpha = .12, position=position_dodge(width=.1)) +
  geom_errorbar(data=total,size=1, aes(x=x0,ymin=CI.lower,ymax=CI.upper,color=type),
                width=0,position=position_dodge(width=.1))+
  geom_point(data=total,aes(x=x0,y=coef,fill=type,color=type), size=2, position=position_dodge(width=.1))  +
  geom_hline(yintercept=0,linetype=2, size=.8) +
  scale_color_manual(values=c("plum4", "black"), name="") +
  scale_fill_manual(values=c("plum4", "black"), name="") + ylim(-.06,.04) +
  ylab("Marginal Effect of ln(Pct.Staff + 1) on Approval")+ xlab("Policy Disagreement") + theme_minimal() +
  theme(legend.position="bottom",
        plot.title = element_text(hjust = 0.5))  + ggtitle("Regular News Followers") 


staff_approval_nonews <- ggplot(total_linear_nonews, aes(X)) +
  geom_line(aes(y = ME, color=factor(type)), size=0.8,alpha = .12, position=position_dodge(width=.1), linetype=1) +
  geom_ribbon(aes(ymin = lower.CI.95.., ymax =upper.CI.95..,fill=type), alpha = .12, position=position_dodge(width=.1)) +
  geom_errorbar(data=total_nonews,size=1, aes(x=x0,ymin=CI.lower,ymax=CI.upper,color=type),
                width=0,position=position_dodge(width=.1))+
  geom_point(data=total_nonews,aes(x=x0,y=coef,fill=type,color=type), size=2, position=position_dodge(width=.1))  +
  geom_hline(yintercept=0,linetype=2, size=.8) +
  scale_color_manual(values=c("plum4", "black"), name="") +
  scale_fill_manual(values=c("plum4", "black"), name="") +ylim(-.06,.04) +
  ylab("Marginal Effect of ln(Pct.Staff + 1) on Approval")+ xlab("Policy Disagreement") + theme_minimal() +
  theme(legend.position="none")+
  theme(axis.text.y = element_blank(),
        axis.title.y= element_blank(),
        plot.title = element_text(hjust = 0.5)) + ggtitle("Non-Regular News Followers") 



staff_approval <- grid_arrange_shared_legend(staff_approval_news, staff_approval_nonews,  ncol=2)
ggsave("FigureH1.pdf", plot=staff_approval,width = 12, height = 5, dpi=600)

###################
# Placebo Test    #
###################

#Run the regression
reg7 <- felm(approval_binary ~ votematchmean2*(traveledonce + travlast + travforward + copartisan + electionyear +
                                                 seniority + nokken_poole_dim1 + chair +
                                                 power + race + female + age)| member_county +
               party_year | 0 | icpsr, data=data,
             weights=data$weight_cumulative)

#Pull out the necessary coefficients to graph
current_interact <- data.frame(coef = reg7$coefficients[20], se = reg7$cse[20],  model="Any Trips This Year x Policy Disagreement")
past_interact <- data.frame(coef = reg7$coefficients[21], se = reg7$cse[21],  model="Any Trips Last Year x Policy Disagreement")
future_interact <- data.frame(coef = reg7$coefficients[22], se = reg7$cse[22],  model="Any Trips Next Year x Policy Disagreement")
current <- data.frame(coef = reg7$coefficients[2], se = reg7$cse[2],  model="Any Trips This Year")
past <- data.frame(coef = reg7$coefficients[3], se = reg7$cse[3],  model="Any Trips Last Year")
future <- data.frame(coef = reg7$coefficients[4], se = reg7$cse[4],  model="Any Trips Next Year")

results <- as.data.frame(rbind(current_interact, past_interact, future_interact, current, past, future))
results$model <- as.factor(results$model)

#Graph it
placebo <- ggplot(data=results,aes(x=model,y=coef,ymin=coef-qnorm(.975)*se,ymax=coef+qnorm(.975)*se))+
  geom_errorbar(width=0,size=1.25,position=position_dodge(width=.8))+
  geom_errorbar(data=results,aes(x=model,ymin=coef-qnorm(.95)*se,ymax=coef+qnorm(.95)*se),
                width=0,size=1.9,position=position_dodge(width=.8))+
  geom_point(size=4,position=position_dodge(width=.8))+
  ylab("Respondent Approval")+
  geom_hline(yintercept=0,size=1.25,colour="indianred",linetype=2) + xlab("") +
  theme_minimal() +guides(fill="none") + theme(legend.title = element_blank()) + coord_flip()

ggsave("FigureE1.pdf",placebo,width = 10, height = 10, dpi=600) 

#####################################
# Regressions Separated by Issue    #
######################################

#County level
reg_civilrights <- felm(approval_binary ~  votematchmean2_civilrights*(log_traveled + copartisan + electionyear + seniority + 
                                                                         nokken_poole_dim1 + chair + power  + race + female + age) | member_county + party_year | 0 | icpsr, 
                        data = data, weights=data$weight_cumulative)

reg_energy <- felm(approval_binary ~  votematchmean2_energy*(log_traveled + copartisan + electionyear + seniority + 
                                                               nokken_poole_dim1 + chair + power  + race + female + age) | member_county + party_year | 0 | icpsr, 
                   data = data, weights=data$weight_cumulative)

reg_foreigntrade <- felm(approval_binary ~  votematchmean2_foreigntrade*(log_traveled + copartisan + electionyear + seniority + 
                                                                           nokken_poole_dim1 + chair + power  + race + female + age) | member_county + party_year | 0 | icpsr, 
                         data = data, weights=data$weight_cumulative)

reg_govopt <- felm(approval_binary ~  votematchmean2_govopt*(log_traveled + copartisan + electionyear + seniority + 
                                                               nokken_poole_dim1 + chair + power  + race + female + age) | member_county + party_year | 0 | icpsr, 
                   data = data, weights=data$weight_cumulative)

reg_health  <- felm(approval_binary ~  votematchmean2_health*(log_traveled + copartisan + electionyear + seniority + 
                                                                nokken_poole_dim1 + chair + power  + race + female + age) | member_county + party_year | 0 | icpsr, 
                    data = data, weights=data$weight_cumulative)

#Note that senator level controls at the congress level are not included for these regressions since they are constant between the two years this issue is asked
reg_crime  <- felm(approval_binary ~  votematchmean2_crime*(log_traveled + copartisan + electionyear + race + female + age) | member_county + party_year | 0 | icpsr, 
                   data = data, weights=data$weight_cumulative)

#Note that senator level controls at the congress level are not included for these regressions since they are constant between the two years this issue is asked
reg_international  <- felm(approval_binary ~  votematchmean2_international*(log_traveled + copartisan + electionyear + race + female + age) | member_county + party_year | 0 | icpsr, 
                           data = data, weights=data$weight_cumulative)

reg_macro  <- felm(approval_binary ~  votematchmean2_macro*(log_traveled + copartisan + electionyear + seniority + 
                                                              nokken_poole_dim1 + chair + power  + race + female + age) | member_county + party_year | 0 | icpsr, 
                   data = data, weights=data$weight_cumulative)

#Save the results
civilrights <- data.frame(coef = reg_civilrights$coefficients[18], se = reg_civilrights$cse[18], model="Civil Rights")
energy <- data.frame(coef = reg_energy$coefficients[18] ,se = reg_energy$cse[18], model="Energy")
foreigntrade <- data.frame(coef = reg_foreigntrade$coefficients[18], se = reg_foreigntrade$cse[18], model="Foreign Trade")
govopt <- data.frame(coef = reg_govopt$coefficients[18], se = reg_govopt$cse[18], model="Government Operations")
health <- data.frame(coef = reg_health$coefficients[18] ,se = reg_health$cse[18], model="Health")
crime <- data.frame(coef = reg_crime$coefficients[14], se = reg_crime$cse[14] ,model="Crime")
international <- data.frame(coef = reg_international$coefficients[14] ,se = reg_international$cse[14], model="International Affairs and Foreign Aid")
macro <- data.frame(coef = reg_macro$coefficients[18], se = reg_macro$cse[18], model="Macroeconomics")

#MSA level
reg_civilrights_msa <- felm(approval_binary ~  votematchmean2_civilrights*(log_traveled_msa + copartisan + electionyear + seniority + 
                                                                             nokken_poole_dim1 + chair + power  + race + female + age) | member_msa + party_year | 0 | icpsr, 
                            data = data, weights=data$weight_cumulative)

reg_energy_msa <- felm(approval_binary ~  votematchmean2_energy*(log_traveled_msa + copartisan + electionyear + seniority + 
                                                                   nokken_poole_dim1 + chair + power  + race + female + age) | member_msa + party_year | 0 | icpsr, 
                       data = data, weights=data$weight_cumulative)

reg_foreigntrade_msa <- felm(approval_binary ~  votematchmean2_foreigntrade*(log_traveled_msa + copartisan + electionyear + seniority + 
                                                                               nokken_poole_dim1 + chair + power  + race + female + age) | member_msa + party_year | 0 | icpsr, 
                             data = data, weights=data$weight_cumulative)

reg_govopt_msa <- felm(approval_binary ~  votematchmean2_govopt*(log_traveled_msa + copartisan + electionyear + seniority + 
                                                                   nokken_poole_dim1 + chair + power  + race + female + age) | member_msa + party_year | 0 | icpsr, 
                       data = data, weights=data$weight_cumulative)

reg_health_msa <- felm(approval_binary ~  votematchmean2_health*(log_traveled_msa + copartisan + electionyear + seniority + 
                                                                   nokken_poole_dim1 + chair + power  + race + female + age) | member_msa + party_year | 0 | icpsr, 
                       data = data, weights=data$weight_cumulative)

#Note that senator level controls at the congress level are not included for these regressions since they are constant between the two years this issue is asked
reg_crime_msa <- felm(approval_binary ~  votematchmean2_crime*(log_traveled_msa + copartisan + electionyear + race + female + age) | member_msa + party_year | 0 | icpsr, 
                      data = data, weights=data$weight_cumulative)

#Note that senator level controls at the congress level are not included for these regressions since they are constant between the two years this issue is asked
reg_international_msa <- felm(approval_binary ~  votematchmean2_international*(log_traveled_msa + copartisan + electionyear + race + female + age) | member_msa + party_year | 0 | icpsr, 
                              data = data, weights=data$weight_cumulative)

reg_macro_msa <- felm(approval_binary ~  votematchmean2_macro*(log_traveled_msa + copartisan + electionyear + seniority + 
                                                                 nokken_poole_dim1 + chair + power  + race + female + age) | member_msa + party_year | 0 | icpsr, 
                      data = data, weights=data$weight_cumulative)


#Save the results
civilrights_msa <- data.frame(coef = reg_civilrights_msa$coefficients[18], se = reg_civilrights_msa$cse[18], model="Civil Rights")
energy_msa <- data.frame(coef = reg_energy_msa$coefficients[18] ,se = reg_energy_msa$cse[18], model="Energy")
foreigntrade_msa <- data.frame(coef = reg_foreigntrade_msa$coefficients[18], se = reg_foreigntrade_msa$cse[18], model="Foreign Trade")
govopt_msa <- data.frame(coef = reg_govopt_msa$coefficients[18], se = reg_govopt_msa$cse[18], model="Government Operations")
health_msa <- data.frame(coef = reg_health_msa$coefficients[18] ,se = reg_health_msa$cse[18], model="Health")
crime_msa <- data.frame(coef = reg_crime_msa$coefficients[14], se = reg_crime_msa$cse[14] ,model="Crime")
international_msa <- data.frame(coef = reg_international_msa$coefficients[14] ,se = reg_international_msa$cse[14], model="International Affairs and Foreign Aid")
macro_msa <- data.frame(coef = reg_macro_msa$coefficients[18], se = reg_macro_msa$cse[18], model="Macroeconomics")

#Combine the results into a dataframe
results<- as.data.frame(rbind(civilrights, energy, foreigntrade,
                              govopt, health, crime, international, macro, civilrights_msa, energy_msa, foreigntrade_msa,
                              govopt_msa, health_msa, crime_msa, international_msa, macro))
sample <- as.data.frame(c(rep("County",8),rep("Metropolitan Statistical Area",8)))
results <- bind_cols(sample, results)
colnames(results)[1] <- "sample"

results$model <- str_wrap(results$model, width = 15)

#Graph it
issuearea_visits <- ggplot(data=results,aes(y=coef,x=model,ymin=coef-qnorm(.975)*se,ymax=coef+qnorm(.975)*se,color=sample))+
  geom_errorbar(width=0,size=1.25,position=position_dodge(width=.7)) +
  geom_errorbar(data=results,aes(x=model,ymin=coef-qnorm(.95)*se,ymax=coef+qnorm(.95)*se,color=sample),
                width=0,size=1.9, position=position_dodge(width=.7))+
  geom_point(size=4, position=position_dodge(width=.7)) +
  ylab("Coefficient Estimate of Policy Disagreement x ln(Visits + 1)")+ xlab("Policy Area") +
  geom_hline(yintercept=0,size=1.25,linetype=2,alpha = .2)+
  scale_color_manual(values=c("plum4", "black"), name="") +
  theme_minimal() +guides(fill="none") + theme(legend.title = element_blank())  +
  theme(legend.position="bottom",
        plot.title = element_text(hjust = 0.5))


ggsave("Figure6.pdf",issuearea_visits,width = 13, height = 6, dpi=600) 

#Now do the same for staff

#County
reg_civilrights <- felm(approval_binary ~  votematchmean2_civilrights*(log_staffed + copartisan + electionyear + seniority + 
                                                                         nokken_poole_dim1 + chair + power  + race + female + age) | member_county + party_year | 0 | icpsr, 
                        data = data, weights=data$weight_cumulative)

reg_energy <- felm(approval_binary ~  votematchmean2_energy*(log_staffed + copartisan + electionyear + seniority + 
                                                               nokken_poole_dim1 + chair + power  + race + female + age) | member_county + party_year | 0 | icpsr, 
                   data = data, weights=data$weight_cumulative)

reg_foreigntrade <- felm(approval_binary ~  votematchmean2_foreigntrade*(log_staffed + copartisan + electionyear + seniority + 
                                                                           nokken_poole_dim1 + chair + power  + race + female + age) | member_county + party_year | 0 | icpsr, 
                         data = data, weights=data$weight_cumulative)

reg_govopt <- felm(approval_binary ~  votematchmean2_govopt*(log_staffed + copartisan + electionyear + seniority + 
                                                               nokken_poole_dim1 + chair + power  + race + female + age) | member_county + party_year | 0 | icpsr, 
                   data = data, weights=data$weight_cumulative)

reg_health  <- felm(approval_binary ~  votematchmean2_health*(log_staffed + copartisan + electionyear + seniority + 
                                                                nokken_poole_dim1 + chair + power  + race + female + age) | member_county + party_year | 0 | icpsr, 
                    data = data, weights=data$weight_cumulative)

#Note that senator level controls at the congress level are not included for these regressions since they are constant between the two years this issue is asked
reg_crime  <- felm(approval_binary ~  votematchmean2_crime*(log_staffed + copartisan + electionyear + race + female + age) | member_county + party_year | 0 | icpsr, 
                   data = data, weights=data$weight_cumulative)

#Note that senator level controls at the congress level are not included for these regressions since they are constant between the two years this issue is asked
reg_international  <- felm(approval_binary ~  votematchmean2_international*(log_staffed + copartisan + electionyear + race + female + age) | member_county + party_year | 0 | icpsr, 
                           data = data, weights=data$weight_cumulative)

reg_macro  <- felm(approval_binary ~  votematchmean2_macro*(log_staffed + copartisan + electionyear + seniority + 
                                                              nokken_poole_dim1 + chair + power  + race + female + age) | member_county + party_year | 0 | icpsr, 
                   data = data, weights=data$weight_cumulative)

#Save the results
civilrights <- data.frame(coef = reg_civilrights$coefficients[18], se = reg_civilrights$cse[18], model="Civil Rights")
energy <- data.frame(coef = reg_energy$coefficients[18] ,se = reg_energy$cse[18], model="Energy")
foreigntrade <- data.frame(coef = reg_foreigntrade$coefficients[18], se = reg_foreigntrade$cse[18], model="Foreign Trade")
govopt <- data.frame(coef = reg_govopt$coefficients[18], se = reg_govopt$cse[18], model="Government Operations")
health <- data.frame(coef = reg_health$coefficients[18] ,se = reg_health$cse[18], model="Health")
crime <- data.frame(coef = reg_crime$coefficients[14], se = reg_crime$cse[14] ,model="Crime")
international <- data.frame(coef = reg_international$coefficients[14] ,se = reg_international$cse[14], model="International Affairs and Foreign Aid")
macro <- data.frame(coef = reg_macro$coefficients[18], se = reg_macro$cse[18], model="Macroeconomics")

#MSA
reg_civilrights_msa <- felm(approval_binary ~  votematchmean2_civilrights*(log_staffed_msa + copartisan + electionyear + seniority + 
                                                                             nokken_poole_dim1 + chair + power  + race + female + age) | member_msa + party_year | 0 | icpsr, 
                            data = data, weights=data$weight_cumulative)

reg_energy_msa <- felm(approval_binary ~  votematchmean2_energy*(log_staffed_msa + copartisan + electionyear + seniority + 
                                                                   nokken_poole_dim1 + chair + power  + race + female + age) | member_msa + party_year | 0 | icpsr, 
                       data = data, weights=data$weight_cumulative)

reg_foreigntrade_msa <- felm(approval_binary ~  votematchmean2_foreigntrade*(log_staffed_msa + copartisan + electionyear + seniority + 
                                                                               nokken_poole_dim1 + chair + power  + race + female + age) | member_msa + party_year | 0 | icpsr, 
                             data = data, weights=data$weight_cumulative)

reg_govopt_msa <- felm(approval_binary ~  votematchmean2_govopt*(log_staffed_msa + copartisan + electionyear + seniority + 
                                                                   nokken_poole_dim1 + chair + power  + race + female + age) | member_msa + party_year | 0 | icpsr, 
                       data = data, weights=data$weight_cumulative)

reg_health_msa <- felm(approval_binary ~  votematchmean2_health*(log_staffed_msa + copartisan + electionyear + seniority + 
                                                                   nokken_poole_dim1 + chair + power  + race + female + age) | member_msa + party_year | 0 | icpsr, 
                       data = data, weights=data$weight_cumulative)

#Note that senator level controls at the congress level are not included for these regressions since they are constant between the two years this issue is asked
reg_crime_msa <- felm(approval_binary ~  votematchmean2_crime*(log_staffed_msa + copartisan + electionyear + race + female + age) | member_msa + party_year | 0 | icpsr, 
                      data = data, weights=data$weight_cumulative)

#Note that senator level controls at the congress level are not included for these regressions since they are constant between the two years this issue is asked
reg_international_msa <- felm(approval_binary ~  votematchmean2_international*(log_staffed_msa + copartisan + electionyear + race + female + age) | member_msa + party_year | 0 | icpsr, 
                              data = data, weights=data$weight_cumulative)

reg_macro_msa <- felm(approval_binary ~  votematchmean2_macro*(log_staffed_msa + copartisan + electionyear + seniority + 
                                                                 nokken_poole_dim1 + chair + power  + race + female + age) | member_msa + party_year | 0 | icpsr, 
                      data = data, weights=data$weight_cumulative)

#Save the results
civilrights_msa <- data.frame(coef = reg_civilrights_msa$coefficients[18], se = reg_civilrights_msa$cse[18], model="Civil Rights")
energy_msa <- data.frame(coef = reg_energy_msa$coefficients[18] ,se = reg_energy_msa$cse[18], model="Energy")
foreigntrade_msa <- data.frame(coef = reg_foreigntrade_msa$coefficients[18], se = reg_foreigntrade_msa$cse[18], model="Foreign Trade")
govopt_msa <- data.frame(coef = reg_govopt_msa$coefficients[18], se = reg_govopt_msa$cse[18], model="Government Operations")
health_msa <- data.frame(coef = reg_health_msa$coefficients[18] ,se = reg_health_msa$cse[18], model="Health")
crime_msa <- data.frame(coef = reg_crime_msa$coefficients[14], se = reg_crime_msa$cse[14] ,model="Crime")
international_msa <- data.frame(coef = reg_international_msa$coefficients[14] ,se = reg_international_msa$cse[14], model="International Affairs and Foreign Aid")
macro_msa <- data.frame(coef = reg_macro_msa$coefficients[18], se = reg_macro_msa$cse[18], model="Macroeconomics")

#Combine the results into a dataframe
results<- as.data.frame(rbind(civilrights, energy, foreigntrade,govopt, health, crime, international, macro,
                              civilrights_msa, energy_msa, foreigntrade_msa,
                              govopt_msa, health_msa, crime_msa, international_msa, macro))
sample <- as.data.frame(c(rep("County",8),rep("Metropolitan Statistical Area",8)))
results <- bind_cols(sample, results)
colnames(results)[1] <- "sample"

results$model <- str_wrap(results$model, width = 15)

#Graph it
issuearea_staff <- ggplot(data=results,aes(y=coef,x=model,ymin=coef-qnorm(.975)*se,ymax=coef+qnorm(.975)*se,color=sample))+
  geom_errorbar(width=0,size=1.25,position=position_dodge(width=.7)) +
  geom_errorbar(data=results,aes(x=model,ymin=coef-qnorm(.95)*se,ymax=coef+qnorm(.95)*se,color=sample),
                width=0,size=1.9, position=position_dodge(width=.7))+
  geom_point(size=4, position=position_dodge(width=.7)) +
  ylab("Coefficient Estimate of Policy Disagreement x ln(Pct. of Staff + 1)")+ xlab("Policy Area") +
  geom_hline(yintercept=0,size=1.25,linetype=2,alpha = .2)+
  scale_color_manual(values=c("plum4", "black"), name="") +
  theme_minimal() +guides(fill="none") + theme(legend.title = element_blank())  +
  theme(legend.position="bottom",
        plot.title = element_text(hjust = 0.5)) 


ggsave("FigureH2.pdf",issuearea_staff,width = 13, height = 6, dpi=600) 


